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NON-ADIABATIC TRANSITIONS THROUGH TILTED AVOIDED 

CROSSINGS 
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(*^ Abstract. We investigate the transition of a quantum wave-packet through a one- 

O^ dimensional avoided crossing of molecular energy levels when the energy levels at the 

crossing point are tilted. Using superadiabatic representations, and an approximation of 
the dynamics near the crossing region, we obtain an explicit formula for the transition 
wavefunction. Our results agree extremely well with high precision ab-initio calculations. 
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r~| 1. Introduction 

-(— > 
C^ The photo-dissociation of diatomic molecules is one of the paradigmatic chemical reactions of 

^ quantum chemistry. The basic mechanism is that a short laser pulse lifts the electronic configura- 

'— ' tion of the molecule into an excited energy state. The nuclei then feel a force due to the changed 

, configuration of the electrons, and start to move according to the classical Born-Oppenheimer 

^ dynamics. Then, at some point in configuration space, the Born-Oppenheimer surfaces of the elec- 

^sD tronic ground state and first excited state come close to each other, leading to a partial breakdown 

^N of the Born-Oppenheimer approximation. As a result, with a certain small probability the electrons 

^^ fall back into the ground state, facilitating the dissociation of the molecule into its atoms. This 

. important mechanism is at the heart of many processes in nature, such as the photo-dissociation 

r^ of ozone, or the reception of light in the retina |19| . For further details on the general mechanism 



o 
o 



X 



we refer to [T^ . 

The mathematical problem associated with photo-dissociation are non-adiabatic transitions at 
avoided crossings in a two-level system, with one effective spatial degree of freedom. Thus, we 
study the system of partial differential equations 

isdtib = Hip, (1.1) 

with V'G L'^{R,C'^), and 

H = -^dll + V{x). 



Above, / is the 2x2 unit matrix, and, with ctx and dz the Pauli matrices as defined in (3.2) 



Vix) = Xix)a^ + Z(x)f7, + dix)I = (^Jg 5ifi)) + di^)I 

is the real-symmetric potential energy matrix in the diabatic representation. Units are such that 
h = 1 and the electron mass TOci = 1. e'^ is the ratio of electron and reduced nuclear mass, 
typically of the order 10~^. The timescale is such that the nuclei (with position coordinate x) 



move a distance of order one within a time of order one. The motivation of ( 1.1 ) and its relevance 
for photodissociation is discussed further in [3]. 
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There is a natural coordinate transformation of ( 1.1 ) that exploits the scale separation provided 
by the small parameter e. The corresponding representation is called the adiabatic representation, 
and is given as follows: Let Uq{x) diagonalize V{x) for each x, and define i!o{x) = {Uq%Ij){x) = 
ip{Uo{x)). Then ipo solves 

ie^tVo = Hoi^o, (1-2) 

with Hq given to leading order by 

fj ^_^g2j^ fpix) + d{x) -eKi{x){ed^)\ ,^ ^. 

° 2 ^ \sKi{x){sdx) -p{x) + d{x))' 

Here, p = \J X"^ + Z^ is half the energy level separation, and 

Z'X - X'Z 



Z2+X2 



is the adiabatic coupling element. A consequence of the choice of time scale in ( 1.2 ) is that solutions 



will oscillate with frequency 1/e. Thus the operator edx is actually of order one. However, we 



have still achieved a decoupling of the two energy levels in ( 1.3 ), up to errors of order e, as long as 
X'^{x) + Z'^{x) > 0. Generically, this inequality is always true: assuming that the entries of V are 
analytic in the nuclear coordinate x, then eigenvalues of V do not cross [T5], and so their difference 
2\/X-^ + Z'^ remains positive. An avoided crossing is a (local or global) minimum of p{x), which 
results in nonadiabatic transitions between the adiabatic energy levels. 

The problem of photodissociation, or more generally of non-radiative decay, can now be formu- 



lated mathematically: assume that (1.2) is solved with an initial wave packet ■i/'in G L'^(R,C'^) that 
is fully in the upper adiabatic level (i.e. the second component of V'in is zero). This is the situation 
just after the laser pulse brings the electrons to their excited state. Assuming that the initial 
momentum is such that the wave packet travels past an avoided crossing, we wish to describe the 
second component of ^o{x,t), to leading order, long after the avoided crossing has been passed. 
By doing this, we predict not only the probability of a molecule dissociating, but also the quantum 
mechanical properties (momentum and position distribution) of the resulting wave packet. 

The difficulty in solving the above problem is that the resulting wave packets are typically 
very small, namely exponentially small in e. As an example, let us assume that the initial wave 
packet "0111 has L^-noim of order one, and that the parameters are such that the L^.^orm of the 
transmitted wave function is expected to be of order 10~^, which we will later see is a fairly typical 
value. This means that any straightforward numerical method with an overall error of more than 
10~^ will produce meaningless results, and thus if we were to apply a standard method (like Strang 



splitting) on the full equation (1.2 1, we would have to use ridiculously small time steps. To make 



things worse, the solution is highly oscillatory. Thus, even though (1.2) is a system of just 1-1-1 



dimensional PDE's, it is not at all trivial to solve numerically. Efficient numerical methods to solve 



(1.2) will therefore require insight into the analytical structure of the equation. 

In [jy, we used superadiabatic representations in order to obtain such insight. We derive a 
closed- form approximation to the transmitted wavefunction at the transition point, which is highly 
accurate for general potential surfaces and initial wavepackets whenever d{x), the trace of the 
potential, is small, but deteriorates when d{x) is moderate or large at the transition point. In 
general, it can not be taken for granted in real world problems that d{x) is small. Therefore, in 
this paper we treat a potential with an arbitrary trace. Our result is weaker than the one in [T]. 
While in the latter paper, we could allow arbitrary incoming wave functions as long as they were 
semiclassical, we essentially require the incoming wave packet to be either Gaussian or a generalized 
Hagedorn wave packet in the present work. However, in that case we still obtain a closed form 
expression for the transmitted wave function at the transition point, and the accuracy is as good 
as in [1,. 
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The importance on nonadiabatic transitions has resulted in much effort to understand them. A 
simplification of the problem is to replace the nuclear degree of freedom by a classical trajectory. 
This approach is both long-established "W and well- understood [11[71[S] and leads to the well-known 
Landau-Zener formula for the transition probability between the electronic levels. This formula 
underpins a range of surface hopping models [l7[ [TOj |9]. Although these and other trajectory- 
based methods ^16J yield reasonably accurate transition probabilities, they are unable to accurately 
predict the shape of the transmitted wave packet 9\. An improvement to the Landau-Zener rates 
is achieved by Zhu-Nakamura theory [12] , which is based on the full quantum scattering theory of 
the problem. However, once again only the transition probabilities are treated, and not the wave 
packet itself. 

It is worth noting that, due to the complexity of the full quantum-mechanical problem of tran- 
sitions at avoided crossings, there are few existing mathematical approaches. The most relevant 
approach to this work is that of [Sj where another formula is given (and proved) for the asymp- 
totic shape of a non-adiabatic wavefunction in the scattering regime at an avoided crossing. For 
simplicity we do not state it here, see Theorem 5.1 of [8j. For comparison, their result looks very 
different to ours, and it is expected that they will not agree in the limit of small e as theirs is 
asymptotically correct, whereas we have aimed for a simple formula which works well for a wide 
range of physically relevant parameters. Nevertheless, our approach is much better suited for prac- 
tical purposes than the formula of [H] , which requires one to calculate complex contour integrals of 
the analytic continuation of some function V that is defined only implicitly. 

2. Computing the non-adiabatic transitions 

In this section we will give a concise overview of our method for computing non-adiabatic 
transition wavef unctions, and explain the various parameters entering the final formula. The 
justification of our method, some extensions and a numerical test will be given in the remainder 
of the paper. 

The data of our problem consists of two parts, the potential energy matrix V and the initial 



wave function. More precisely, we assume that we are given p{x) and d{x) as in (1.3 1, and that 
p has a unique global minimum in the region of space that we are interested in. We choose the 
coordinate system such that this minimum occurs at x = 0, and we thus have 

p{x)^d + 0{x^), d{x) = do + Xx + 0{x^). 

The transmitted wavefunction only depends on on A and p, but unfortunately the latter quantity 
does not enter in a simple way. Under the reasonable assumption that the matrix elements X and 
Z are analytic functions of a; at least close to the real axis, then so is p^. We write p{q)^ — 5'^+g{qY 
where g is analytic and g{Q) = 0. Since g^ is quadratic at 0, a Stokes line (a curve with Im(/9) = 0) 
crosses the real axis perpendicularly, and, for small 5, extends into the complex plane to two 
complex zeros of p, namely qs and g|. We define, for any complex z, the 'natural scale' [1] 



r(z) = 2 / p(OdC, 

and write t^ = T{qs), where qs by convention is the complex zero with positive imaginary part. 
We write 

Tr = Re (ts), Tc = Im (ts), 

which are the two parameters that enter into the transition formula. When we are given p in a 
functional form, neither the computation of its complex zeroes not of the complex line integral 
leading to Tg is a problem numerically, and can be carried out to any required accuracy. How- 
ever, in the case of radiationless transitions, the potential energy surfaces are often known only 
approximately. As our final formula will depend very sensitively on the value of rg, small errors 
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in this quantity will lead to wrong predictions. This is not a fault of our method, but a general 
obstruction to any numerical method aiming to calculate small nonadiabatic transitions: namely, 
since our formula below agrees very accurately with ab-initio computations, and depends very 
sensitively on t^, getting p wrong will lead to wrong results regardless of the method used. In a 
way, it should not be too surprising that when looking for a very small effect, we need to get the 
data right with very high accuracy. But it does pose a serious practical challenge when trying to 
predict small nonadiabatic transitions. 

As for the initial wavefunction, first of all we assume it to be initially concentrated in the upper 



electronic energy band. This means that we will consider equation (1.2 1 with initial condition 
'00(2;, 0) = ((/)+(a:, 0),0)"^. The restriction of this work, when compared to the case with A = 
considered in [1] , is on the form of 0+ (a; , 0) , which we require to be either Gaussian, or a finite linear 
combination of Gaussians, or a Hagedorn wavefunction. For the present exposition, we restrict to 
the case where it is Gaussian. The first step of our algorithm is straighforward: 

Step 1: Solve the upper band adiabatic equation \edt4'+ = H^ijj^^ ■0+(O) = 0+('iO), where 
if + = —£^9^/2 + p{x) + d{x). This can be done either by direct Strang splitting, or using the 
theory of Hagedorn wave packets [TT]. For a transition to occur, we need the wave packet to cross 
the transition region near a; = 0, where p is minimal. So we monitor the expected position {X) = 
J x\'il>^{x)\'^ dx and stop the evolution when (X) = 0, say at time ig- Let us write (j){x) — ip+{x, Iq). 
4> is Gaussian up to errors of order e [TT], and centered at x = 0. Thus we have 

?^(A;) = exp(-^(fc-po)') (2.1) 

with parameters po (the mean momentum) and c. Above, we used the semiclassical Fourier trans- 



form, cf. (3.8) 



Given these initial data, we need to define one further derived quantity. Put tlq — -^, where 
ko is part of the solution of the pair of equations 

k = ^rf + 45, r; = fc(l - l£too) ) . (2.2) 

Again, the numerical value of uq is easy to obtain. In what follows, we will use the abbreviation 

rC ^Tf{k) = v/fc2 - 4(5. 



Step 2 Put 



0_ (fc) w — , ^^exp 



a2,oao 1 + "o,2ai o ~ "i,oao,i"i,i 



af 1 - 4a2,oao,2 



(2.3) 



>4(5, 



with 

sgn(fc)rc + irr - rfriQe ~ 4:cS{r]* - po) 



"1,0 



2Sy/e k + Tj* 



2(no + l)e'/^X . , 2(no + 1)A£ 

"O'l "" IT, — * ' ai.i = -ir^ +-77— — ;^T^> (2-4) 

k + 7]* (k + rj*)'' 

2Snoe + r]*^ e . 25A 2(no + 1)A2£ 

"2.0 = w^ C-— -— — — -, ao9 ^ ■ 



and 



8(52 2(fc + r7*)2' ""'^ (k + ri*) [k + i]*)^ 

(no + l)^eAao(5 1 / uqS \ ,. ^'^ 

^(^") = ~2(nn + l)2A2s2 + 2(52a2 " 2 "'^^*"" I (n^^TD^J + '^"^^^"^ 4 
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where above uq — \/pq + 45 +po- While formula (2.3 1 is trival to implement on a computer and 



produces accurate results, it would of course be desirable to interpret the various terms in a phys- 
ically meaningful way. However, we have been unable to do this. On the other hand, everything 
except the factor 6""^*^^") jg obtained by approximations and exact Gaussian computations. The 
latter factor is more tricky. We include it because we found a discrepancy in the phase of the 
transmitted wave function between our formula without the factor and numerical ab-initio calcu- 
lations, which is probably due to one of our approximations below being too crude. This phase 
discrepancy is removed by essentially computing the phase in the case of an incoming Gaussian 
when the parameter c diverges, i.e. infinitely small momentum uncertainty, which gives ip{po), and 
subtracting that. While this fixes the discrepancy with the numerics, we do not, as yet, fully 
understand why it does so, and where the original inaccurate approximation has been made. 

One important property of fipo) is that it is constant in both k and x and hence will not 
affect any quantum mechanical expectation values. It would however play a role when we consider 
interferences. 

The final step of our algorithm is again straightforward: 
Step 3 Solve the lower band adiabatic equation with initial condition 0_ , i.e. solve 

where H^ = —e'^d'^/2 — p{x) + d{x), and 0_ is the inverse semiclassical Fourier transform of 
01. For times so large that V'- has support far away from the transition region, it describes the 
transmitted wave function of equation (1.2 1 with great accuracy, see section ro] 



3. Evolution in the Superadiabatic Representations 

3.1. Superadiabatic Representations. The key idea for deriving our transition formulae is 
to study the evolution in a suitable superadiabatic representation. For a careful discussion of 
the theory of those representation, we refer to [3]. Here we give only some intuition and the 
mathematical facts. The n-th superadiabatic representation is implemented by a unitary operator 



Un acting on L^(M, C^), and its main property is that it diagonalizes the right hand side of (1.1) 
up to errors of order £"+^. Thus, the adiabatic representation (1.3) is the zeroth superadiabatic 
representation, and in general 

^\^2r , (P{x)+d{x) £"+iX+ 



where K^ are the n-th superadiabatic coupling elements. They are usually pseudo-differential 
operators, and so are the C/„. The useful consequence of switching to the superadiabatic represen- 
tation is that now the evolution of the second component "0" of V'n — Un'ip, subject to ip' (~oo) = 0, 
is given by 

^-(i) = -ie"r e-i(*-^)^"i^„-+ie-i^^^0ds, (3.1) 



up to relative errors of order e. Thus, provided we can control K~ , (3.1) gives the transmitted 
wave function in the n-th superadiabatic representation to high precision. 

There are some apparent problems with this idea. Firstly, it is far from clear how we hope to 
control K~ . Secondly, the superadiabatic unitaries are in general very hard to calcualte, and as 
such this formulation does not allow the adiabatic wavefunction to be easily obtained. Thirdly, we 
have to decide which value of n we want to use. The sequence K~ is expected to be asymptotic 
in n, so after initially decaying rapidly (in an appropriate sense) it will start to grow beyond all 
limits when n is taken to infinity. The second problem is resolved when we study the wavefunction 
in the scattering regime, well away from the avoided crossing. In this case, for potentials which 
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are approximately constant, it is known that C/„ and Uq agree up to small errors depending on the 



derivatives of the potential |15) . and (3.1 1 can be used to calculate the transmitted wavefunction. 
For the value of n, in [3] we showed for a special choice of parameters p, k that there exists an 
'optimal' n for which ipni^) builds up monotonically, corresponding to a single transition. This 



n is given by the set of nonlinear equations (2.2 1 that we have seen in the previous section. We 



expect this set of equations to hold in general, and have obtained very good results by using it 
here. 

The problem of calculating K^ turns out to be reducible to a set of differential recursions, which 
we will now give. The discussion follows the one in [7 very closely, the only difference being that 
we now include a nonzero trace d{x) in the Hamiltonian. All the calculations and arguments are 
almost the same as in [3], so we will omit them. 

We change from the spatial representation to the symbolic representation (see e.g. [TS]) by 
replacing a; by g € M and ied^ by an independent variable p £ R, where the factor e takes into 
account the semiclassical scaling. We need to introduce some further notation: we rewrite the 
potential as 



t/('?)-P(9)(^gi^(0(^)) cos(%))j+°'^'^^1^0 1 




which defines 9{q). It follows that the unitary transformation to the adiabatic representation is 
given by 

Uo{q) = 

Hence the Pauli matrices in the adiabatic representation are given by 

crx(g) = Uo{q)a^Uo{q), (7y{q) = Uo{q)ayUo{q), CT^{q) = Uo{q)aJJo{q), 

where 

^0 l\ /O -i\ A 



1 oy "y"Vi oy ""~vo -1/' ^^-^^ 

and we have used that U^ — Uq. 

A direct calculation confirms, with 1 the 2x2 identity matrix: 

Lemma 3.1. We have 

d'^Viq) = a,,{q)(^M + KiqWM + c„(g)l, 
where an{q), bn{q) and Cn{q) are given by the recursions 

ao{q) = p{q), a«+i(g) = a'^,{q) + 0' {q)bn{q) 

boiq) = 0, fe„+i(<7) = b[^iq)-d'{q)a^{q) (3.3) 

coiq) = d{q) Cn+i{q) = c'„{q) 

We then have the following explicit recursion for the coupling elements: 
Theorem 3.2. The Hamiltonian in the n-th superadiabatic representation is given by 

H^{e,p,q) 2 1+ l^£"+i«;^^(p,5) -p{q) + d{q) ) + \p{e^^+^) 0{e') 



where 



K^+i(p,g) = ~2p{q){xn+i{p,q) ±yn+i{p,q))- 
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Setting X nip, q) — J2m=oP"^ "^^^nil) ™^^ similar expressions for yn, Zn andwn, the coefficients 
x™ to ui™ are determined by the following recursive algebraic- differential equations: 



. ^'(9) 1 



X™ = zT = < = 0, m = 0, 1, y1 = -i^' yi = (3-4) 

with 

( Lm/2J 



/or n o(i(i, and 
7;™ — - 

i/n+l ^ 




J(«)'-0'zr)-2 Y. ^[r.^^-r^^,^(^-a,y::;^r+h,u,::;^r+c,x::-^^^) 

7 = 1 ^ ' 



, L™/2J 



(2i)J 



Lm/2J 



-('7/;"M'-2 Y^ -^ /n+l-m+jw m-2j , , m-2i , m-2ix 

a^n J ^ 2^ (21)-?'^ J' j>>'^J^n+l-j +Oj^«+l-j +CjUI„^;^_jJ, 



for n even. The coefficients a„ to Cn are given by Lemma 3.1 



Proof. The proof is analogous to those of Theorem 3.4 and Proposition 3.5, with some easy alter- 
ations due to the presence of d{x) . D 

We note that, as in the trace- free case in [3], y™ = for all in when n is even and x™ = 
z™ = w™ — for all m when n is odd. Furthermore, from the above equations, it is obvious that 
a;;j* = y™ = z™ = w™ = for odd to. 

We now have an explicit expression for k,^ and may therefore also calculate K~ , the superadi- 
abatic coupling element, which is the Weyl quantization of the symbol k~: 

K^^ix) ^^ f d^dy«± (^,0 ei«--^)^(y), 

and from the recursions in Theorem 13.21 it is clear that 



where the k„^„_j can be calculated explicitly. Determining the asymptotics of this two-parameter 
recursion is a very tricky problem to which we have no solution. However, in the regime of large 
p (meaning, large incoming momentum) the sum is well- approximated by the j = n term. For 
p — 0{e^^^^), this can be made rigorous on the level of the superadiabatic Hamiltonian, while a 
full asymptotic investigation of the transitions in this regime is still work in progress [5]. Here, we 
use this approximation without further justification, and find that it gives good results even for 
relatively small values of p. 

The asymptotics of the term k~q can be determined explicitly in the following generic case. 
Without loss of generality, we assume that the avoided crossing occurs at a; = 0, specify the initial 
wave packet at i = 0, and write p{q)^ = (5^ + g{q)^ where g is analytic and ^(O) = 0. As is 
standard in asymptotic analysis (see e.g. (4j), the asymptotic behaviour of k~q is determined by 
the complex zeros of p. Since g^ is quadratic at 0, a Stokes line (a curve with Im(p) = 0) crosses the 
real axis perpendicularly, and, for small 5, extends into the complex plane to two complex zeros of 
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yO, namely qs and qg. As argued by Berry and Lim, in the natural scale T{q) = 2 J„ p(q)dq and near 

q = 0, the adiabatic coupling function has the form Ki{q) = ^^^ ( m - * TT^ ^ i^riTio))), 

with Ts = T{qs). In particular, k^ has no singularities fro \t\ < jr^l, and no singularities of order 



^ 1 for \t\ ^ jr^j. As can be seen from Theorem 3.2 solving the recursions for k„ requires taking 



high derivatives of ki. By the Darboux principle, the asymptotics are dominated by the complex 
singularities closest to the real axis, Tg and t|. Hence, to leading order, we find 

-,7,0(9) = ^P(g)(n- l)!((^4^ - (^3^). (3.5) 

Using the definition of the Weyl quantisation, a direct calculation [3 shows 



^n^o -EL- {§y{dinl,ix))i~ied,r-=. (3.6) 

3.2. Approximation of the Adiabatic Propagators. In order to determine a closed form 



approximation for (3.1), it is necessary to approximate the adiabatic propagators. This is in 
contrast to the situation in [3j where the model was chosen such that p is constant, and thus the 
adiabatic evolutions were trivial in Fourier space. 



The first insight is that the operator K~q given in (3.6) is sharply localized: K'^f will only be 
significantly different from zero if either / or some of its derivatives have some support overlap with 
K^O' which means they must be concentrated near the real solution of Re (T(g)) = Re {ts) that is 
closest to g = 0. We will refer to this solution as the transition point. In Section [4?T] we will see that 
relevant values of n are of the order 1/e; furthermore, for large n we have (l + x^)"" « e~"^ , and 
so K~Q and its derivatives are concentrated in a i/e neigbourhood of the transition point. Since the 
time scale is chosen such that the semiclassical wave packets (which have width of order y/e) travel 
at speed of order one, the dominant transitions come from a time interval of order ^/e around the 
transition time, which we define to be the time when the expected position of the incoming wave 
packet crosses the transition point. 

Let us pick a coordinate system so that that the transition time is s = 0. We cannot, however, 
choose the transition point to be at x = 0, since we have already fixed x = to be the local minimum 
of p. On the other hand, one of our later calculations relies on the fact that the transition point 



is at least in a V£ neighbourhood of 0, see Section 3.3 So from now on, we will always assume 
that the transition point does indeed have this property. This assumption can be justified by the 
observation that for sensible potentials, the real and imaginary parts of the complex zeroes of p are 
coupled, and are either both relatively small or both large. However, in the latter case, transitions 
tend to be so small that they are physically uninteresting. That said, it would of course be much 
preferable to be able to treat arbitrary transitions, but we cannot do this yet. In what follows, 
we will always pretend that the transition point is a; == 0, although for the calculation in the next 
paragraph below this is not yet strictly necessary. 

The above considerations allow us to replace the potential in the full adiabatic dynamics by its 
first Taylor approximation, as the following formal calculation shows. We take Hf^ :— —e^d'^/2 ± 

6 + Xx and wish to show that e e — e e ^ is small. We have 

--sH^ --sH7= --sHf ( -sHf --sH^ ^\ --sHf [ c, ( -sHf --sH^\-, 

ee — ee i=ee Mee ^ee — lj = ee ^ Orie^ ^ee )dr 

zsH^(i,u± TJ±\\ „-hsH^ 







We now note that H^ — H^ is quadratic near a; = and hence the integrand is of order 1 in a y/e 
neigbourhood of zero. Hence the left hand side is bounded by the length of the integration region 
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y^ and, to leading order, it sufBces to replace (3.1 1 by 



(3.7) 



where we have not altered the s-independent propagator. 

We now find it convenient to switch to the Fourier representation by applying the semiclassical 
Fourier transform 

f{k) = ^ /e-i'=V(g)dg=^/(f). (3.8) 



We define Kn through Kni'^ = Knip , and a direct calculation ^3| gives 

I -^e 

Fourier transforming both sides of (3.7), we see that i\jn is given by a double integral: 



Ipn {k,t) 



: e e ^ ' 



i^Txe 



As 



dr;ei^^^"W<+i,o(^-'7)(^) 



n+l i 



-sffiWIe 



(^), 



where H^ (H^) are the approximate (exact) adiabatic propagators in momentum space. 
By the Avron-Herbst formula, the approximate propagators are given exactly by 



In particular, we have 



isHtik) 



zsH^ik) 



^sH+iv) 



= e se e 



Xsdk 



^((k^±2S)s-\ks^) 



e^ e-^^^'^ ei(^''-'*)^ e^^^'' 



(3.9) 



iA^s^ 



e 6e e " e 26 



Xsdj^ 



— {t]^-2S)s p^Aijs 



e2e' 



In order to make use of these expressions w e m ust understand the effects of the shift operators, 
where e'^*^'' f{k) = f{k + Xs). Using (3.9) in (3.7) we note that, due to the invariance of the integral 

under j] t-^ rj — As, we may apply the rj shift to the left with opposite sign. Hence k^_|_]^ q (fc — r/) 
is unaffected and {k + 77)"+^ 1-^ {k + i] - 2As)"+^ 

Shifting the remaining propagator in k by —As, the remaining multiplicative parts of the prop- 
agators are given by 

exp [^([(fc - As)2 - 2d]s + \{k - As)s2 - {rj^ + 26)s + Xt]s^)] . 

Simplifying this expression and inserting it into (|3.7|) gives 



Tpn (fc,t) 



le" 



/27re 



^^i^C^) / ds dr^ik + r,- 2As)"+i«-+i^o (fc - 77) 



{{k^-rt^-iS)s-{k-ri)Xs^) ^e ^^y 



X e2e 



(3.10) 



3.3. Fourier transform of the cou plin g el eme nts. In order to make use of (3.10), we require 
the Fourier transform of n^o- Using (^3.8| on (3.5) gives 

\n+l 



^nfi \k) 



1 



/27r£ 

1 
/27re 



e-i'^''J^^p(q)(n-l)! 

IT 
i i^+l 



1 



1 



27r 






dq 



{t-t;Y {t-ts)' 



where we have used dr — 2p{q)dq. 
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It is now that we need the transition point to be at or near q = 0. Provided this is so, we can 
use that p has a minimum (5 at g = 0, and expand q{T) = ^ + 0{t^). Note that no second order 
term is present. As the remainder of the integrand is concentrated in a ^/e neighbourhood around 
g = 0, we keep only the first order term, giving 



'^n.O W 



1 



/27re 



e 2Se 



kr 1 



n+1 



27r 



■{n-iy. 



1 



1 



i^-^s)" {r-TsY" 



dr. 



T^ 



i-i 1 



and hence 



We now note that , ^ ,„ = (-1)""^7 

1 ;n+l r j 

^"-1 ' e"23l'=^a"-i 



/2^ 27r 



(-1)' 



-2ir^ 



(0 



T - TS 



dr. 



dr 



Using the identities pik) = ^/(f), dff{k) = (ifc)"/>), /(a;-a)(fc) 



-ia/c 



f{k) and the 



standard Fourier transform 



4^W 



-a|fc| 



gives 



'*n,0 (^) ~ 1 



.72(5 1 /A:\"-i _^ 

^(2(5)" vi; 



e 25el*''l e '2Se'' , 



where we have used ts = Tr + iTc- Inserting this formulation into (3.10[) gives 



V'n ik,t) 



--tH-(k) 



ine 



e e 



d. /dr,(fc + ^)(l-|A|)"+i(iL_^) 



X e 25, 



Ss!''-''! e"53i('^-'') g2i(('='-'''-4'5)«-('=-';)^'*') 0e(^). (3.11) 



4. Evaluation of the integral 



4.1. The choice of n. Equation (3.11 ) still depends on the parameter n, the order of the supera- 



diabatic representation. For choosing n, we attempt to use the same argument that was employed 
in [3] in order to obtain universal transition histories. The idea then and now is that the modulus 
of the integrand in (3.111 depends on n, while the phase does not. We will thus try to choose n 



such that stationary phase and maximal modulus occur at the same point, making it possible to 
perform asymptotic analysis on the integral. 

We recall the assumption that Tr is small, and consider the imaginary part of the exponent. 
Indeed, we will set Tj. = in what follows. This simplifies the analysis and does not seem to greatly 



affect the accuracy of the final result. Differentiating the phase of (3.11 ) with respect to s and r] 
gives 

(4.1) 

(4.2) 



(fc^ - 77^ _ 4j) _ 2A(fc - 77)s = 
-27/s - Xs^ = 0. 



Note that if A = then there is only one solution, namely k^ — rj^ = AS and s = 0; this remains a 
solution if A 7^ 0. 



For a simultaneous solution to (4.1 ) and (4.2 ) (i.e. stationary phase for both integrals) we require 
either s — Q and k^—rf'—AS = 0, or As = — 2?7 and rj the solution to —brf+Akrj+k^ — AS = 0. In the 
second case, for k = 0(1), we see that 77 and hence s are also of order 1. We have already discussed 
that we expect the significant transitions to occur only when s — ©(e^/^), and we therefore expect 
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this solution to contribute only a negligible amount to the transmitted wave packet. So from the 
stationary phase condition, we obtain s = and k^ — rf —45 — 0. 



For the modulus, we assume the case of a Gaussian wave packet of the form (2.1 1. Differentiating 



the logarithm of the modulus with respect to ry and s and equating to zero leads to the equations 

2A 



(,1 + 1)- 



M'n-vo)- i^v 



ne 



277 



fc2 



•q^ 



- {n + l)£ 



V- 

2Xs 



k-2i 



(k + T]){k + T] - 2\s) 



= 0, 



= 0. 



(4.3) 
(4.4) 



Equations (4.1 )-(^4.4| cannot be solved simultaneously, which shows an interesting difference of the 



present case when compared to the non-tilted case treated in [T] and [3]. To make progress, we 
argue that the choice of the optimal superadiabatic representation should depend only weakly on 
the trace A of the potential. Therefore, we allow A to vary as well as n, 77 and s, and obtain the 
joint solution s = A = 0, and n and 77 fulfilling n 
future always use this value of n, denoted tiq. 



eko 



with ko the solution of (2.2 1. We will in 



4.2. Rescaling. Recall that the wavepacket moves a distance of order 1 in time of order 1, and, 
for a semiclassical wavepacket, is of width of order e^/^. Hence for times of order e'' with 7 < 1/2, 
in position space, the wavefunction is localised well away from the transition region. It follows 
that there should be little contribution to the integral outside s e [~e'^,£'^] for 7 < 1/2. We thus 
restrict the s-i ntegra l to this region. 

We rewrite (3.11) as j^ exp{—^tH^ (k)) J^dij J_^-, ds g{k,r],s) with 



5(77, fc, s) = exp 



nlo. 



'tl-rf 



2Xs 

k+7) 



2i5£l 



log(fc + 77) + (71 + 1) log I 
^'Mk -v) + ii [{k' - rf - AS)s - Xik - r^)s'\ cjy^{r,). 



We now note that, in order for the phase of the integrand to be stationary in s, we expect rj « 
77* — zt\/k^ — 45. For a semiclassical wave packet which has sufficient momentum to move past 
the avoided crossing, the choice of sign will correspond to the sign of the mean momentum of cf)'^ . 
For this choice of 77* to make sense, it is clear that we require k^ — 45 ^ 0, and so introduce the 
cutoff function Xk^>4S- The physical meaning of this cutoff is clear when one considers 77 to be 
the incoming momentum and k the outgoing momentum: since the potential gap is 2(5, by energy 
conservation we have k'^ /2 = rf' /2 + 25, and since we require 77^ > for the wave packet to move 
past the crossing we have fc^ > 4(5. 

We now set rj = fje^'^ + 77*, where 77 is of order 1 and rescale the s integral by s = se^'^, which 
causes the domain of the s integral to be at least of order 1, and tend to the whole of M as e — >■ 0. 



Using /jj d77 /f ^^ ds g{k, t],s) ^ e /jj d77 /f 
from now on, we are interested in 



-1/2 



.1/2^ 



girls'^' +vi*,k,e^'^s) 
+ {n + 1) log (1 - 



exp 



77, log (1 



1/2 ds5(A:,e^/^77 + ?;*,£ 



45 



2\se 



1/2 



25c 



k+ri'+rie^'^ 

+ h [(-^^e - 27777*ei/2)^gi/2 _ ^(^ 
We now discuss the evaluation of these two integrals 



+ log(fc 

1/2 



^/^s), and removing the tildes 

1/2) 



— 7/ — rje 

7/* — ?7e"^' ^)s^e] 



'2* 



■q 
■{k 



{e^l\ 



rj* — rj£ 
+ 77 



1/2N 



(4.5) 



4.3. The s integral. Since the wave function (j)' is independent of s, we now aim to perform the 
s-integration explicitly. We now consider the regime where e is small and k is of order 1. This is 
necessary as we wish to expand the logarithm term in powers of s, and require that , ^ /_f — j^j <IC 1. 
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This holds since, from the hmits of integration, we see that at worst se^'"^ ^ e* with 7 > and 
77* ^ fc ^ 77 ^ 1 . Expanding to second order gives 



log 1 



2Ase 



1/2 



1\SE 



1/2 



2A^s^e 



fc+77*+r;ei/2/ 



fc+))*+i)ei/2 (fc+))*+i(ei/2)2 



(4.6) 



In the small-e limit, e^ ^ " — >■ c», which, combined with the above expansion reduces the s-integral 
to a Gaussian integral of the form 

/ exp(as^ + /3s)ds = J exp ( ) , for Re (a) < 0. 

Jr \ a \ 4a/ 



In this case, we have 
(where Re (a) < 0) and 



(fc+r,*+ei/2^)2 2\'^ n } ^ 2 '/ 



ii 



2A(n+l)£^/'' _ • * _ m!^ 2 
"fc + ^.+£l/2^ l?? '/ 2 '^ ■ 



It therefore remains to calculate the integral over r/ 

4A/7re 



iei/^A 



-1/2 



X exp 
X exp 



log (1 - 2!^±?mV^) + iog(fc + ,7* + ,^£1/2) _ ^|fc _ ^* _ ^,1/2 1 - i^(fc - rf - r;£V2)- 

/ 2A(n+l)5V^ . . ie^^%2\V 2(n+l)A^s jA /,. *^ , ie^^^A \ 



For a general 0^, we can say little else, and the integral must be computed numerically. However, 
in the important case where c^f is a Gaussian, we can derive a closed-form approximation, which 
is in excellent agreement with the full dynamics. The main idea is to approximate the integrand 



in (3.11 1 in such a way as to produce a Gaussian integral. The first hinderance to this comes from 



the log terms, which we now consider. 

4.4. Expansion of log terms. Along with the expansion in ( |4.6[ ), we have 

log(fc + ry* + ,7£i/2) = log(fc + rf) + log (l + ^) « log(fe + ,7*) + 2£^ - .^ig^, 

log (1 - ^z!£±2m:£!^) « _:Pe±2gr£^ „ _^(^4^2 ^ 4^3^*^3/2 ^ 4^2^*2^)^ 

where we have once again used that k,jf ,1] ^ 1. 

In order to produce a Gaussian integral, it is necessary to make a number of justifiable approx- 



1/2„1-P_ 



P 



1,2 in (4.6) around rye^/^ = and neglecting terms 



imations. Expanding (fc + 77* + e ' rj) 

of order larger than e in all three logarithm expansions reduces them to: 



log 1 



2Asei/^ 



2As£i/2 



2As£i) 



2A^s^e 



fc+r;*+r;ei/2; ~ fe+r,' T^(fe+^.)2 (fc+,,.)2! 

log(fc + rj*+ ^eV2) « log(fc + rf ) + 2£^ - 5^^, 

^"& l-^ 45 J ~ 45 8.52 ■ 



(4.7) 



Note that all three expansions now contain terms of at most order two in s and 77 and thus are of 
the form required for a Gaussian integral. 
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4.5. Explicit closed form. One final simplification is necessary to obtain a Gaussian integral 

i^rfs and ^^ 



(4.5) still contains the third order terms, namely ^-^^^r] s and „'i/2 ?7s ■ But here we recall that 



the staionary phase argument required s = 0, and in the scaled variables also 77 = 0. This allows 
us to remove the above terms: not only are these terms already the highest order in e, but since 
we expect the main contribution to the integral to come from the region around (s, 77) = (0, 0), the 
effects of these terms is negligible. 



Inserting the expansions (4.7) into (4.5), ignoring the third order terms in s and 77, and setting 
0^(77) to be the Gaussian 0'^(ry) = exp(— 1(?7 — po)^) gives, for 77 sufficiently small, 

g{k,rie'^/^ +r]*,e^/^s) = exp(a2,o?y^ + ai.o?? + "i.i^s + "o.is + ao,2S^), 



th the Qfij given in (2.4). Note that the sgn(A:) in ai.o is necessary if we wish to deal with 

77 > 0. 



negative momenta: for fc > 0, we have k — rj* > and hence, for small e, k — if 



^1/2 



Therefore \k 



V 



.1/2 



ri\ 



1/2 



|fc — 77* — e^/^771 ~ \k — ■i]*\ + e^i^rj. 
Gaussian integration now gives 



\k — ri*\— e^^'^rj. However, for fc < we have k — rj* — e^^^rj < and 



drjds g{k,rie ' +?7*,e ' s 



.l/2_«^ ^ 



27r 



y 4a2,oao,2 - a?,i 



exp 



"2, otto 1 + ao,2ai - ai,oQ:o,iQ^i,i 



Off 1 - 4Q;2,oao,2 



(4.8) 



which holds for Re ( -^^ 

\ Q2,0 



4ao,2) > 0. 

2 
We now check that the above constraint Re (,^^ — 4ao,2) > is satisfied for a suitable parameter 

regime. For ease of analysis, we note that no is approximately given by Tc/(e\/pQ + 46) = 0{s^^). 
Taking e to be small, to leading order we find 

"2,0 — Afi aX2 L . 



TcV 



ai 1 



-ir; 



2nn£A 



(fe+r,')2 



-177 



; 25\ 

ao.2 = -i(7T+;^ 



2(no + l)A^£ 
(fe+rr)2 



I 2tcX 

{k+rr)^y/pl+45 

^ ; 2<5A 2r<;A^ 



Note that the real part of — 4ao,2 is non-negative, so we need only check the sign of Re {af i/a2.o)- 
Using al ^ — —rj* 



iri\' 



4rcAi)* 



Re(^-4ao2)> 



gives 



7f^{k+-n')\pl+iS)-irlX' 



^Q2.o -">-/' {k+ri*)i^pl+iS L S5-^c^pl+iS+2ST,+T,pl 

Since Tc > 0, this is clearly positive when po is sufficiently large. Hence the regime of interest is e 



small and po large. We then have 



V'« {k,t) 



'-tH- 



1 



2w4a2,oQ;o,2 - a?^i 
X {rj* + k)e 



: exp 



"2,0^0.1 + 0^0,2^1.0 - "1,0^0,10^1,1 



in'^poV e-2fel''-"*l e 



'25. 



"1,1 ^ 4a2,oao,2 



with the aij as given in (2.4 1 



We note that setting A = gives ao,i 
form (see [3]) 



Xk^>4S, (4.9) 

c^o.2 — and ai^i — i?]*, and returns the n-independent 



Tpn (fc,0) 



iv'+k) -liv'-pof e-^|fe-r,*| e-'sfcC^-"*) 

2|J7'| 



Xfe2>4(5- 
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4.6. Asymptotics of aij. We now show that, under suitable assumptions, the ai_j given in (2.4) 
may be somewhat simphfied. Note that uq — 0{e^^) and so Q2,o == ^(1) ^-iid ai.o — 0{e^^^. 



However, the two terms which come from the n-independent prefactor k + rj in (3.11 1 are of lower 
order than the remaining terms in the respective a^.j, and hence, for small e may safely be neglected. 
From the point of view of exponential asymptotics this is completely natural; one would normally 
fix the slowly varying terms (i.e. independent of e, and in this case of n) at the stationary value 
of the integrand. For clarity, we now have 

"2,0— 45 852 L, ai,0 — 25 eVs "T 25ei/2 ^" ^SeVa ■ 

One additional simplification is possible when pq is large and the wavefunction is quickly decay- 
ing (i.e. c is also large). In this case the modulus of the integrand is negligible unless rj* is close 
to Pq. For such a range, using ng ~ Tc/ {e y/p^+iS) , shows that the first three terms in ai^o above 
are all negligible. Further, if the potential is symmetric, t^ = and we may use the approximation 
Q^i.o — 0. In addition, in this limit, and with the assumption that A is not too large, we see that 



the second terms in each of ai^i and ao,2 in (2.4) are negligible. To conclude, for e small, po and 
c large and A not too large, we have 

"2.0- S5^^-'' "^'°"2Vl' "M«-i^, 

2A i25X 



^{k + r)r,*' "'' (k + r,*)- 



4.7. Additional Phase shift. While testing the formula ( |4.9[ ) against ab-initio numerics, we 
found a discrepancy by a phase shift which, in the region where the wavefunction has significant 
magnitude, is constant in k. We believe that this effect comes from one of the approximations 
detailed above, but have currently been unable to determine its exact cause. For many applications 
this phase shift is unimportant. Since it is constant in fc, all expected values of observables are 



correctly reproduced by (4.9 1 in the case of a single Gaussian wave packet. Where the phase shift 
begins to matter is for interference phenomena, and when considering a superposition of Gaussians 
(see below) such that their centres are at significantly different locations in k; then the phase shift 
will not be constant in k any more, and we will get wrong predictions for position expected values. 
It is therefore desirable to have a method of removing this effect of the approximations. We 
now describe a heuristic method which has proven to be effective for a wide range of potentials 



and initial Gaussian wave packets. Consider (4.9) for A 7^ and the wavefunction normalized by a 



prefactor ■y/c/(7re). Note that if A = the following argument is invalid. However, setting A = in 



(4.9) we see that the phase depends only on r^, which agrees with that of IJ and the corresponding 
numerics. 

We are going to consider the phase of the transmitted wave function in the limit c — > 00, 
i.e. the incoming wave packet approximating a ^-function at rj = po- Since the numerical phase 
shift is independent of k, we need to choose a value of k at which to evaluate this phase. In the 
classical picture, from energy conservation we see that the the transmitted wave packet should be 
approximately a 6- function at fc = yPg ~^ '^'^' ^^^^ hence we consider this value of k, where the 
sign of the square root is chosen to match that of po- 

We are therefore interested in 



'^n (y; 



p2 + 4,5,0)^ 


. ^ 1 


"a2,oao,i + "o,2"i,o " ai.oao.iaia" 


^ — cxp 

V^e 2^4a2.oao,2 - a?,i 


L "1,1 ~ 4a2,oao,2 J 




xaoe2fcl'-"l e2Se^«Xk^>4S, 
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where aq = \l v\ + 4(5 + po and 60 = y v\ + 4(5 — po- We now investigate the phase of this wave 
packet when c ^- 00 and note that there are contributions from both the square root and the 
exponent. 

We write a2,o = /32,o — c, and, since r\* — po, this is the only term that depends on c. Consider 
first the prefactor: 



Vc 



I 



y 4Q;2,oao,2 - af ,1 y 4(/32,o - c)(3;o,2 - a?,i y-4<3!o,2 + ^(4^2,o"o,2 - a?,i) 



\/-4ao,; 



(4.10) 



For the exponent, we have 



0:2,00:0,1 + oo,2a?,o ~ oi,oao,iOi,i 
"1 1 ~ 4(3!2,oao,2 



il32,0 - c)(3!q ;^ + (3^0,20? ~ 01,000,101,1 



*1.1 



4(/32,o - c)ao,: 



^0,1 



^(/32,OOo 1 + 0:0,20? - Oi,oOo,lOi,i) 



a, 



0,1 



4oo,2 + c("l 1 ~ 4/32,000,2) 



4a, 



0,2 



(4.11) 



It remains to determine the phases of (4.10) and (4.11). We write ai^o = /3i,0j Oi^i = /3i.i+i7i,i, 
oo,i = /3o,i and ao,2 = /?o,2 + 170,2, with Aj,7i,j e M. For ( |4.10[) , we note that -q;o,2 = -(;5o,2 + 
170.2) =: r e'^ , where 9 — arctan ( ^^^ J , giving the phase of (4.10 ) as — ^ arctan ( 3^ j . Using that 



70,2 = -A&0/2, /3o,2 = -2(no + l)X'^e/al andao^o = 45 , this is -^ arctan yj:^;;^:p[j^ 



For the (4.11) we have 



QQI? 



l)£^/'^X/ao, this phase is given by ~ 2in+iy\^e^+25^a'„ 
(no + l)^£Aao(5 1 



4(;3o,2+i7o,2) iil3l2+7i,2) 



Hence, using that /3o,i — 2(no + 



2(no + l)2A2e2 + 2(52a2 2 



arctan I 



aoS 



{?2o + l)eA 



25e 



bo. 



One further adjustment seems to be necessary. One would expect that the phase is continuous 
in A, and we know that for A = 0, the phase is —^bo. However, the limit of the A-dependent 
terms in above expression is — l/2arctan(sgn(A) sgn(ao)oo) = — sgn(A) sgn(po)7r/4 and hence we 
take the phase shift to be 



(no + l)'^e\ao6 1 

^(^") = -2{no + irX^s^ + 2S^al ' 2 "'■^*"" 



apS 
(no + l)eA 



sgn{Xpo)' 



(4.12) 



which seems to give very good numerical results for a wide range of all parameters. 

To summarize, we now have an explicit closed form for the transmitted wave packet given an 
initial Gaussian of the form (2.1 1: 

- — ^ e 



ztH' 



1 



: exp 



2i/4(3!2,oOo,2 - ai,i 



02,00o 1 + Oo,20i - Oi^oOo,lOi,i 



a 



1,1 



402,oOo,2 



X (7;* +fc)e~M^I'^-"*l e-'^"'^''*) e-''^(P°)(^^(7y*)xfc2 



>4(5, 



(4.13) 



with ^p{po) as given in ( |4.12[ ) and the aij as in (2.4), or, alternatively, with the simplifications 
from Section 4.6 ng is given as indicated in (2.2). We thus have finished the justification of the 
algorithm given in Section [21 
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4.8. Phase shift for large momentum. For large momentum, we use the approximations qq = 
2po and uq = Tc/{epo) and get the phase shift (f{po) as (p{po) « ~^ rg?2+4?^pg ~ ^ arctan ( ^j + 
sgn(Apo)f Note that if po ^ oo then (p(po) -^ 0. More concretely, we are interested in the rate at 
which it goes to zero when we write Pq in terms of e. Letting pq = e~°' gives 

1 T^XS 



If: 



e T^A^ei 



4^2gl-3a 



1 Me-^"5 
— arctan ^ — 

2 V TrX 



sgn(Apo)- 



Hence if a > 1/3 we see that ip{e^°') — >■ as e — ?► 0. We note that this value of 1/3 is the same 
value as that for which we have rigorous bounds on the errors [2]. 

From this analysis, it appears that the phase shift is a consequence of taking momenta that are 
too small (or equivalently, e that are too large). 



5. Non-Gaussian incoming wavefunctions 

5.1. Extension to Hagedorn v^ravefunctions. We note that a general Hagedorn wavefunction 
[Sj is a Hermite polynomial multiplied by a Gaussian. By linearity of the integral, it is sufficient 
to consider the case (p'^irj) = tjP cxp (~c{rj ~ PqY /e), p £ N. We perform the same rescaling as in 
Section 4.2 and note that the monomial prefactor becomes {rj£^/'^+'q*Y — X]?=o (^)(?7e^^^)-'?7*^^~-'^ 
Using t 



le same arguments as above, we obtain for each j the integral 

dr/ds (e^/^Ty)-' exp(a2,o?7^ + ai,o?y + ai,i»7s + ao,iS + ao,2S^) 



We now note that 9^ exp(ai.o»7) = ?]■' exp(ai_o'7) and since differentiation with respect to ai^o 
commutes with the integral, we have 

dryds {e^^^rjY exp(a2,o'7^ + aifiV + ai^irys + ao,is + ao,2S^) 



- ■ 27r 




"1,0 / 

w4a2,oao,2 - 


-<i 


1 2^ 


= exp 



exp 



a2,oao,i + "o,2ai,o - ai,oao,iai,i 



4a2,oao,: 



4a2,oao,2 



"2,oao 1 



4Q2,oao,2 



di 



ai.o exp 



ao,2ai,o - ai,oao,i"i,i 



4a2nao.2 



In more generality, we wish to compute d^f where / — exp{—aa^+ba) — exp(— ^(a— ^)^ + ^). 
It is clear that this will be / multiplied by a scaled and shifted Hermite polynomial. In fact, we 
have d^f = {—^/2ayHj{^/2a{a— j^))f, where Hj is the probabilist's Hermite polynomial of order 
j (namely chosen such that the coefficient of the leading order is 1). 

In our case, we have a = —-:^ /^ — and o = ' ■ mT^nrr jl — "'^ ^^ 



-4q2,o"o,: 



2 — -. — ■ , givmg ^ — '"'' ^''' . Hence 

afj^-4a2,oao,2 ' ° <= 2a 2qo,2 



(3.11 ) with 4''^{ri) = rf exp {—c{rj — pq^ /e) is given by 



fpn {k,t) 



\tH-{k) 



Xk^>iS 



'2\/'ia2fiao^2 - aj i 



: exp 



"2, otto 1 + "0,20:1.0 



ai,oao,iai,i 



Q^i.i ~ 4Q;2,oao,2 



;{n'-po) 



X {t^* +k)e-^^'' -P")' e-sftl'^-"*! e-'^siC'-'n') 



P , 
3=0 



s}j^<p-j) 



(«?.! 



2ao,; 



4a2,oao,; 



j72 



X iJ, 



2ao,; 



a^ 1 - 4a2,oao,2 



1/2 



(ai,o 



ao,iai,i 
2ao,2 



with the aij as given in (2.4 1. 
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Using the identity Hp{x 

2qo,2 \1/2/ 



VLo(^)x-^i/,(y), withx 



2ao 



-1/2 



y=(a 



-4q;2.0Q!0,2 ' 



"1,0 



ao,iai.i ' 
2ao,2 ■ 



gives 



-4q2,0"0,2' 



and 



i^ (fc,t)«e-i*^"W 



Xk'^yU 



X (?7* 



X H, 



2W4a2, 0^0,2 
2ao,2 



: exp 



e 25e 



Ifc-';*! „-i 



a2,oaoa + «o,2Q;i,o - ai.oao.iQ^i,! 
Q^i.i ~ 4a2,oao,2 

2ao,2 



25? ('=-''') eP 



(«?i 



p/2 



4a2, 0^0,2 



-1/2 



2q;o.2 



4a2,oao,2 



4a2,oao,2 
1/2 / ao itti.i 



("1:0 



2ao,: 



We are interested in the leading order behaviour with respect to e. From ( |2.4| and using 

Tip = ©(e^^) we see that a2,Ojai,i)Cto,2 are all 0(1) whilst ai^ and ao,i are 0(e~^/^). Hence 

2 _?°'^ = C'(l)i and aifl — "°'^"^''" = 0(6"^/^), which in particular shows that the prefactor 

is 0{e^) whist the argument of the Hermite polynomial is 0{e^^). Thus, to leading order, only 
the highest power of the Hermite polynomial contributes, giving 



i^n {k,t) 



-^tH' 



1 



: exp 



a2,oao 1 + ao,2ai o ~ ai,oao,iai,i 



2w4a2,oao,2 - "la 



4a2 ,0^0,2 



i2i;(fc-')-)^*P^,2>45, 



which is precisely (4.91 with the Gaussian replaced by if exp {—c{i^ — poY /e) 



We note that the error in this closed form is expected to be of order ^/£. Whilst this could 
be improved by taking further terms in the expansion, in the following we choose to concentrate 
on the case of a wave packet which has been decomposed into a linear combination of complex 
Gaussians. The main reason for this is the heuristic phase correction which is discussed in Section 
4.7 From numerical studies, we see that this works well only for Gaussian wave packets, and 
without this correction, the relative error between the formula and the 'exact' numerical wave 
packet is of the order of 10%, compared to an error of around 2% for a Gaussian wave packet with 
the phase correction. 

5.2. General wave packets as superpositions of Gaussians. Due to the strong reliance on 



the wave packet being a Gaussian in the preceding discussion, formula (4.9) is not immediately 



applicable to general wave packets. However, we propose a simple algorithm which allows use 
of (4.9 1. For a given semi-classical wave packet specified on the upper level well away from the 



crossing region, we evolve it using the BO dynamics on the upper level until the mean position of 
the wave packet coincides with the crossing point (which we choose without loss of generality to be 
at a; = 0). We then transform into Fourier space and decompose into complex Gaussians, giving a 
wave packet of the form 

N 



Va, 



/ ^ 1 j exp I 



{P-P^Y 



exp ( 1 



pXj' 



(5.1) 



where in position space Xj is the offset from the crossing point. 

We now need to deal with the fact that the Gaussians reach the crossing point at different times. 
We fist note that, for small e, this should be a small effect for semiclassical wavepackets: since the 
wave packet is localised in a ^/e neighbourhood of zero in position space, we have Xj = 0{-\/£). 



As discussed in Section 3.2 on a y^ neighbourhood of the origin and for times of order -y/e, the 
dynamics are well-approximated by the explicit propagators (3.9). Since the wave packets move 
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with speed of order one, this is stiU the region of interest and we may simply insert the complex 



Gaussian into (3.111 



Applying the rescaling as described in Section 4.2 gives an extra term in the exponent in (4.5) 
of the form ixj{r]* + e^''^r])/e. The rj* term combines with the Gaussian term in ry* to give the 
wave packet evalua ted at 77* as before. The remaining term provides a contribution of the form 
iXj/e^/^ to ai,o in (2.4). 



It is now easy to see that in the small e limit this term is negligible. Since Xj = 0{y/e) we see 
that the new term in ai is order one. In contrast, the dominant terms in ai.o are of order e"^'^ 



and one may apply (4.9) directly to the complex Gaussian. 

We note that for the values of e under consideration in the numerics, ignoring this correction 
increases the relative error by the order of 0.1%, which is quite significant given the high accuracy 
of the final formula. In the implementation of the non-Gaussian wave packet below, we therefore 
included the additional term ixj/e^''^ in the expression for ao,i. 

The above analysis suggests a simple and efficient algorithm for calculating the form of the 
transmitted wave packet, even if not of Gaussian form, given an initial wave packet V'-oo located 
well away from the transition point in position space. 

(1) Evolve the initial wave packet on the upper BO level using the uncoupled BO dynamics 
until its centre of mass reaches the transition point. This can either be pre-determined 
by finding the point at which the two energy levels are closest, or, as would be required 
in higher dimensions, by monitoring the energy gap at the centre of mass over time and 
determine its minimum. 

(2) Transform the resulting wave packet into momentum space and decompose into a linear 



combination of complex Gaussians as in (5.1). 



(3) Apply formula (2.3) to each complex Gaussian in turn and take the corresponding linear 
combination. 

(4) Evolve the resulting transmitted wave packet using the BO dynamics on the lower level, 
until the centre of mass reaches the scattering region. 

Assuming that the energy levels become constant in the scattering regime, the computed wave 
packet will agree up to small errors with that computed using the full coupled dynamics. 

Note that step (2) may be accomplished using a standard numerical recipe such as non-linear 
least squares optimization. In practice the formula is more accurate for narrow Gaussians (since 
this improves a number of the approximations including the choice of fixed ng and the heuristic 
phase shift) and thus it may be worth constraining the variances of the Gaussians. Since the 
application of the formula is cheap (simply multiplications in Fourier space over a region in which 
the modulus of the wave packet is significant - comparable to one time step in uncoupled B-0 
dynamics) and step (3) scales linearly with the number of Gaussians, increasing the number of 
Gaussians whilst decreasing their variances would be a reasonable approach to increase accuracy. 

It is important to realise that, although this algorithm performs a molecular dynamics calcu- 
lation using Gaussian wave packets, it does not share the obstructions of most Gaussian-based 
methods (see e.g. HH). These occur mainly due to the Gaussians being not orthogonal, and the 
resulting ill-conditioning of various matrices under time evolution. Since we only require that the 
wave packet is decomposed into Gaussians at the crossing point, transmitted, and re-summed on 
the lower level, we do not encounter such problems. In fact, one is free to choose any method 
of propagation on the adiabatic levels, for example the method of Hagedorn wave packets [TT] . 
something which will be important in higher dimensions, where simple grid-based methods are 
prohibitively expensive. 



BORN-OPPENHEIMER TRANSITIONS 



19 




Figure 1. The two adiabatic energy surfaces V± = ±p(a;) + d{x) with p{x) = 
y/X'^{x) + Z'^{x). Z = atanh(a:) + (3x^ / cosh{x), X ^S, d^^ Atanh(a;), with the 
parameters a — 0.5, f3 — —0.4, S — 0.5 and A = 1. Note that the avoided crossing 
(minimum of the energy gap) is at a; = 0. 



6. Numerics 



We now compare the results of formula (2.3) to those of high-precision fully-coupled numerics. 
For ease of demonstration, we set the transition point to be at x = and t = and choose to 
specify the initial wave packet 4> as a linear combination of complex Gaussians in momentum space 
at the crossing time. This simplifies the implementation of the above algorithm, as 0^ is already in 
the required form. Further, we set e = 1/50, which gives reasonably small transition probabilities, 
whilst still enabling the 'exact' calculations to be performed. Note that, when transformed to 
position space, both examples have mean position zero. 

To begin both the full numerics and the implementation of the above algorithm, we evolve (f) on 
the upper BO surface to large negative time (i.e. to a position where the potentials are essentially 
flat) to give a good approximation (/)_t ~ i'-oo- 

The full numerics were performed using a symmetric Strang splitting in matlab with initial 
condition (J)-t, which is run to a time i* > where once again the potentials are essentially flat. 
In particular, for times i > i*, the lower component ||'!/'o'(^)ll i^ constant. We then evolve '0o'(^*) 
backwards in time to i = and compare its Fourier transform to formula (2.3 1. The calculation 



was performed on a grid with 16,384 points in both the position ([—40,40]) and corresponding 
momentum ([—12.87,12.87]) spaces, with T — t^ — A and 1000 time steps. Doubling both the 
number of space and time gridpoints produces a wave function which differs from this computation 
by around 0.01% in the L^ norm, and hence we take the numerical simulation to be 'exact'. 

We choose Z = atanh(a;) -I- /3a;^/cosh(x), X = S and d = Atanh(a;). For these choices, 6 
and A correspond to their earlier use, the ratio a'^ /S determines the second derivative of p at the 
transition point, and /3 primarily affects the asymmetry of the potential. In particular, /3 = gives 
Tr = 0. We set a — 0.5, /3 = —0.4, S = 0.5 and A = 1. This leads to the two potential surfaces 
given in Figure [T] with ts = —0.16611 + 0. 537721, which can be easily calculated numerically. 
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Figure 2. Top: The absolute value of the transmitted wave packet ip~ as given 
by (|2.3|) (solid, left axis) and the error as compared to the numerical solution 

<j)~ computed as described using the fully coupled dynamics (dashed, right axis). 
Inset is the initial Gaussian wave packet in momentum space at the transition 
time. Bottom: As in the top plot but showing the argument (phase) of the wave 
packets. 



The first wave packet we treat is given by the complex Gaussian Aexp [—c{p — po)^ /{2e)), 
with po=5, c = 1/(2(7^), a = \/2 and A chosen such that the wave packet is normalized in L^. 



The second case we consider is a linear combination of three complex Gaussians of the form (5.1 1 
where \Aj\ = A, j = 1,2,3, which in turn is chosen to normalize the wave packet. The remaining 
parameters are given, with c = l/(2cr^) by 
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Figure 3. As Figure [2] but for the non-Gaussian wave packet described in the text. 
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In both cases, the relative error is less than 2% over the full interval where the transmitted 
wave function is essentially supported. The transition probability ||V'~IP in both cases is of the 
order 10^^ (3.03 x 10~^ and 3.48 x 10^^ for the Gaussian and non-Gaussian cases respectively). 
In addition to these two examples, we have tested a wide range of parameters for the both the 
potentials and semi-classical wave functions, and all results are good to within a few percent. They 
deteriorate only when e (and thus also ||^~||) becomes too large and we leave the adiabatic regime, 
or when po (and thus also ||V'~||) g^ts too small and our many approximations requiring that po 
is suitably large break down. In particular, the relative error is less than a few percent when the 
transition probability is in the range 10^^-10^^^. 
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